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The observation of a granular structure in high-energy heavy-ion collisions can be used as a 
signature for the quark-gluon plasma phase transition, if the phase transition is first order in nature. 
We propose methods to detect a granular structure by the single-event intensity interferometry. 
We find that the correlation function from a chaotic source of granular droplets exhibits large 
, fluctuations, with maxima and minima at relative momenta which depend on the relative coordinates 

. of the droplet centers. The presence of this type of maxima and minima of a single-event correlation 

' function at many relative momenta is a signature for a granular structure and a first-order QCD 

\ phase transition. We further observe that the Fourier transform of the correlation function of a 

■ granular structure exhibits maxima at the relative spatial coordinates of the droplet centers, which 

can provide another signature of the granular structure. 
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PACS numbers: 25.75.-q, 25.75.Nq, 25.75.Gz 



I. INTRODUCTION AND SUMMARY 

> . 

■ Recently, much progress has been made in the experimental search for the quark-gluon plasma • The occurrence 
of jet-quenching and collective flow suggests the presence of a very dense matter produced in high-energy heavy-ion 
collisions P, 01 ■ A very important question is whether the produced dense matter is a quark-gluon plasma. If it is a 

■ quark-gluon plasma, it will undergo a phase transition from the quark-gluon plasma phase to the hadronic phase. It 
is of great interest to search for the signature for the phase transition of the quark-gluon plasma. 

The signature for the phase transition depends sensitively on the order of the transition. Previously, Witten and 
many other workers noted that a granular structure of droplets occurs in a first-order QCD phase transition, and 
the observation of the granular structure can be used as a signature for a first-order QCD phase transition |^-[l5j|. 
In a recent spinodal analysis, Randrup found that the spinodal instability has a high degree of amplification and 
the instability of the most rapidly amplified wavelengths grow predominantly, leading to spinodal (granular) density 
patterns that may indeed be used as a diagnostic tool for a first-order QCD phase transition ■ Assuming that the 
granular particles are arranged in regular rapidity intervals in momentum space, Randrup studied methods to identify 
the associated momentum clumping using N-particle momentum correlations [l5j . 

We would like to study the granular source in configuration space. Among the many different ways to study 

■ a quark-gluon plasma, intensity interferometry (HBT interferometry) is the best experimental tool to examine the 
\ space-time density distribution of the produced matter It can be utilized to study the granular structure 

that occurs in a first-order phase transition of the plasma. 

It should be pointed out that without a careful study of the phase-transition dynamics and post-transition evolu- 
tion, it is not known at present how much the granular density pattern of the phase transition will remain and become 
detectable by HBT interferometry. It also remains a subject of current research how early a post-transition configura- 
tion the HBT interferometry really detects, as it was shown by quantum treatments of the multiple scattering process 
and the collective flow that HBT interferometry measures the density distribution at a configuration earlier than the 
thermal freeze-out configuration [s^ l40l l4ll |42| . To continue our progress in the search for the phase transition of 
a quark-gluon plasma, it is reasonable to start with the working hypothesis that the density fluctuations that occur 
during a first-order phase transition are so large that some remnants of the granular droplet distribution remain after 
the post-transition evolution, and these remnants form an initial chaotic source of granular droplets. The assumed 
chaoticity of the particle source then leads to the result that the correlation function of two identical bosons in HBT 
interferometry contains information on the space-time density of the source. A granular density distribution of the 
emitting chaotic source then distinguishes itself from other density distributions and should reveal its characteristics 
in HBT interferometry. 

Suggestions to examine the granular structure in connection with the phase transition of the quark-gluon plasma 
have been presented previously 0-01 ■ Recent high-energy heavy-ion measurements give a ratio of iiout/-Rsidc ~ 1 
[45I li^ . which is contrary to most theoretical expectations H^-jl^. A granular emitting source of droplets was 
proposed to explain this HBT interferometry puzzle T^. The suggestion was based on the observation that in the 
hydrodynamical model the particle emission time scales with the radius of the droplet. Particles will be emitted 
earlier if the radius of the droplet is smaller, as in a source of many droplets. An earlier emission time will lead to a 
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smaller extracted HBT radius -Rout- On the other hand, -Rsidc depends on the distribution of droplet centers and is 
relatively independent of the droplet size. It increases with an increase in the width of the droplet center distribution 
and the collective-expansion velocity of the droplets. As a result, the value of i?out can lie close to i?side for a granular 
quark-gluon plasma source. A direct investigation on the granular density structure of the emitting source is, however, 
needed to confirm the occurrence of a granular structure in high-energy heavy-ion collisions. 

Previously, Pratt et at. ^ studied the granular density structure in HBT interferometry by considering pairs of 
identical bosons in a distributed source in which the centers of the droplets are distributed according to a spherical 
Gaussian distribution with a standard deviation Ra/\/2. They calculated the correlation function C(|q|) for a pair 
of identical bosons of relative momentum q from such a distributed source with _Ro = 4 fm. They found that there 
are differences in the magnitudes of the correlation function C(|q|) at |q| greater than about 80 MeV/c in a granular 
structure. They also found that the correlation functions are relatively smooth functions of |q| even for small numbers 
of droplets (2 and 4) . 

The source considered by Pratt et al. in Ref. corresponds to a source in which the centers of the droplets have 
been distributed first and identical bosons are then subsequently examined from such a distributed source. In actual 
dynamics with droplet formation in a first-order phase transition |^ , one envisages the formation of localized droplets 
in each single event following the phase transition, and particles are emitted and evolve from these localized granular 
droplets. Although there can be a distribution of the centers of the localized droplets over many different heavy-ion 
collision events (of similar other global characteristics) , the centers of the droplets in each single collision event can 
be localized. 

Noting that the centers of the droplets in each single event can be localized in space and time, we would like 
to consider a different way to study the granular structure in HBT interferometry. We propose experiments to 
look at HBT interferometry for each event individually, to study its associated density distribution. There are many 
advantages in the single-event HBT interferometry, if it can be carried out with sufficient accuracy. First, the granular 
structure is the result of a large fluctuation of the density distribution in space and time as it undergoes a first-order 
phase transition. The large fluctuation of the density in space and time is encoded into its corresponding Fourier 
transform, whose absolute square gives the correlation function of identical bosons. By considering the correlation 
function of all pairs of identical bosons from this single event, one can obtain pertinent information on the density 
distribution of the event. One can even invert the correlation function with an inverse Fourier transform to obtain 
an integral equation for the source density distribution of the event. Secondly, as a result of the large fluctuation 
in the dynamics of a first-order phase-transition, the granular structure of one event will be quite different from the 
granular structure of another event, even though many other global characteristics can be very similar. By examining 
the space-time structure of each event individually one can study the large fluctuation of different events in the 
phase transition, which is another characteristic of a first-order phase transition. In other words, there are both 
large fluctuations in a single event as well as large fluctuations among different events in a first-order QCD phase 
transition. Finally, any averaging over a set of collision events to obtain an average correlation function, as in the 
use of identical boson pairs from different collision events in a 'multi-event' analysis or in considering bosons from a 
'distributed source', will wash out the large fluctuations in the single-event correlation function. As a consequence, 
the large fluctuations of the correlation function that are inherent in a first-order phase transition in single events do 
not become prominent. A single-event HBT analysis brings out the prominent features of large fluctuations of the 
correlation function. 

At intermediate energies and the low-energy end of high-energy heavy- ion collisions, it is difficult to use the single- 
event HBT interferometry as the number of detected bosons is small in a single event. It becomes necessary to pick 
identical boson pairs from a large number of events with similar global characteristics in a 'multi-event' analysis in 
order to provide sufficient statistics so as to extract useful, but average, space-time characteristics concerning this 
group of collisions. As a large number of events are included in the sampling, the fiuctuation of the correlation function 
that may be present in the single-event HBT interferometry is suppressed by the averaging procedure. 

In a nearly head-on collision at very high energies, the number of identical pions is of the order of a few thousand. 
The number of observed identical pions n,r is a small fraction of this number. For example, the number of identical 
pions detected in the STAR Collaboration in the most central Au-Au collisions at RHIC is of the order of a few 
hundred [s^ . Although the number of pairs of identical pions in the event varies as njrl'T'ir — l)/2, only a small 
fraction of these pairs have relative momentum small enough to be useful in a HBT analysis. Clearly, whether or not 
a single-event HBT measurement can be carried out remains to be tested. If the statistics of identical bosons turns 
out to be insufficient for such an analysis using the present detectors in heavy- ion collisions at RHIC, there remains 
the possibility of performing a single-event HBT analysis at RHIC with detector upgrades or with heavy-ion collisions 
at LHC. It will also be of great interest to carry out a few-event HBT analysis in future work both theoretically 
and experimentally to see how the degrees of fluctuation changes as the number of events increases. The rate of the 
change of the fluctuations will provide information on the underlying fluctuation in single-event and event-to-event 
fluctuations in HBT interferometry. The few-event analysis will require a good theoretical understanding of the 



3 



single-event HBT interferometry. 

We shall show that the correlation function of a granular density distribution has large fluctuations, with maxima 
and minima at locations which depend on the relative coordinates of the droplet centers. These local maxima and 
minima arise from the constructive and destructive interference of identical bosons originating from two different 
droplets. Their interference pattern therefore depends on the coordinates of the droplet centers. The occurrence of 
these maxima and minima in single-event correlation functions is a good signature for the granular structure and a 
first-order QCD phase transition. 

It is desirable to take advantage of the relationship between the density distribution and the correlation function to 
invert the correlation function by Fourier transform to obtain an integral equation for the density function. Previously, 
using the Koonin-Pratt formalism ,19.. .27i|. Brown, Danielewicz, and their collaborators studied methods to invert 
the correlation function by expanding the two-particle source function in spherical harmonics [ssl Wh . li^ . For an 
irregular density distribution as one encounters in granular droplets, an expansion in terms of spherical harmonics 
is inadequate. We have developed a general three-dimensional Fast Fourier Transform method to invert the correla- 
tion function which will facilitate its application to HBT interferometry, should accurate experimental single-event 
correlation functions become available. 

Using the property that the Fourier transform of the correlation function is the folding of the density with itself, 
we shall prove analytically that the Fourier transform of the correlation function of a granular density distribution 
has maxima at spatial locations governed by the spatial relative coordinates of the droplet centers. The presence of 
this type of maxima in the Fourier transform of the single-event correlation function provides another signature for 
the granular structure of the density distribution. We shall demonstrate the presence of these maxima numerically by 
inverting the correlation function with the three-dimensional Fast Fourier Transform (FFT) method. A comparison 
of these results with the exact analytical solution shows the feasibility and the high degrees of accuracy of our 
three-dimensional FFT method. 

This paper is organized as follows. In Section II, we briefly summarize the basic relation between the density 
distribution and the correlation function. We give in Section III the correlation function for a simple periodic structure 
and show that it possesses many correlation function maxima. In Section IV, we study the correlation function of 
a general granular density distribution. We flnd that the correlation function exhibits large fluctuations. We show 
that the locations of the maxima and minima of the correlation function are determined by the relative coordinates 
of the droplet centers. In Section V, we show that the Fourier transform of the correlation function has maxima at 
locations which also depend on the relative coordinates of the droplets centers. The feasibility of the Fourier transform 
method is demonstrated by numerical examples. Section VI concludes our discussions. To facilitate the use of the 
Fourier transform method using experimental single-event correlation functions, (or functions that fit the experimental 
correlation function), the method to invert a correlation function using the 3-dimensional Fourier transform is given 
in Appendix A. The method to obtain the iterative solution of the density distribution from the Fourier transform of 
the correlation function is given in Appendix B. 



In a HBT interferometry measurement, one measures the correlation function P(fci,/c2), the probability for the 
detection of one boson of one 4- momentum fci = (/cio,ki) in coincidence with the detection of an identical boson 
of another 4-momentum /c2- An alternative representation of the correlation is given in terms of the correlation 
function C{ki,k2) defined as the ratio P{ki,k2)/P{ki)P{k2), where P{ki) is the probability of detecting a boson of 
4-momentum fc^. 

Because of the symmetrization of the wave function of identical particles, there is an interference of the two 
histories for two identical bosons to propagate from two source points to two detection points. For a chaotic source, 
this interference leads to a relation between the two-boson correlation function and the Fourier transform of the source 



II. DENSITY DISTRIBUTION AND CORRELATION FUNCTION 



density p{x; fci, fc2) by p^-[4^ 




(1) 



where 




\/.f{ki,x)f{k2,x) 



(2) 




/(fc, x) is the Wigner function, and x = {t, r). We shall assume that the source is indeed chaotic so that Eq. |^ holds. 
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The density function p depends on the average momentum of the pair of particles P = {ki + fc2)/2 and the relative 
momentum q = ki — k2- We shall neglect final-state interactions and assume that the density function p is independent 
of q. Then the average momentum P becomes just a label for the correlation function and the density function, 

= 1 + I /da;e^«Xx,P)|'. (3) 



For brevity of notation, we shall leave the label P implicit so that the correlation function C{q) and the density 
function p{x) in fact refers to those with a pair momentum P. The label P can be re-introduced when it is needed. 
In this implicit notation, the correlation function C{q) is related to the source density function p{x) by 



C{q) = l + \j dxe"'-^p{x)\ . (4) 
We normalized the density function according to 

dx p{x) = 1. (5) 



It is easy to show that the correlation function C{q) is a real function greater than unity and possesses an inversion 
symmetry, 

C{q) = C(-g). (6) 

Because of the simple relationship between C{q) and p{x) in Eq. Q, there is a one-to-one mapping between the 
density function p{x) and the correlation function C{q). The correlation function C{q) will exhibit special properties 
associated with the density distribution p{x). Conversely, information concerning p{x), including its fine structure, is 
encoded in C{q) so that direct information on p{x) can, in principle, be obtained by decoding C{q) (see discussions 
in Section V). 

In single-event HBT interferometry experiments, one measures P{ki,k2) in the same event and the correlation 
function C(fci,fc2) is obtained by dividing the joint probability Piki^k^) by P{ki) and P{k2). Ideally, in the single 
event, if the number of detected particles is large enough, the single-particle probability distribution P{k) can be well 
determined using detected particles in the same event. However, when one is faced with low statistics in the total 
numbers of detected particles, it may be necessary to obtain P{k) from a larger set of events (with similar global 
characteristics) to determine the correlation function C{ki,k2). 



III. AN EXAMPLE OF A REGULAR GRANULAR SOURCE 



It is instructive to study many simple examples to find out the characteristics of the correlation function of granular 
sources in HBT interferometry. We shall begin by considering a regular granular structure. While we do not expect the 
QCD phase transition to lead to a completely regular granular structure, the remnants of the density amplifications of 
the most rapidly amplified wavelengths that grow predominantly during the first-order phase transition may remain, 
and these amplified wavelengths may show up in the underlying gross structure of an otherwise irregular granular 
pattern [T^ . It is thus useful to study the correlation function of a regular periodic granular structure. The simplicity of 
the solution also provides interesting insight into the nature of the correlation function of granular density distributions. 

Accordingly, we consider a spatially periodic granular source of the form 

3 

p{t, r) = fit) YlicosKjX, + 1)9{R - r), (7) 

3=0 

where the granular structure is characterized by an overall spatial dimension R, and the wave numbers Kj (or 
wavelengths dj = 2ti /Kj) with j = 1,2,3. The function f{t) contains a multiplicative normalization constant chosen 
to satisfy the normalization condition (|SJ. The source is a periodic granular distribution contained in a sphere of 
radius R. The Fourier transform of p[x) is 



p(9o,q) - /(9o) / dr e-^i "-a(r)6(r). 



(8) 
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where 



filo) = J 



dt e'^°'f{t), 



a(r) = JJ(cosii'ja;j + 1), 

6(r) = e{R - r). 
Using the folding theorem of Fourier transforms, we obtain 



/ 



dr e-''»""a(r)6(r) 



(03«(q')Kq-q'), 



where {a(q),6(q)} are the Fourier transforms of {a(r),6(r)} respectively. They are given explicitly by 

3 

aid) = {2nf H [6iqj) + 5{qj - Kj)/2 + 5{qj + Kj)/2] , 



6(q) = AttR 



qR 



From these results, we obtain 



P{q) = f{qo) E c(Ai)c(A2)c(A3 



|q - QA1A2A3I-R 



where 



=(A) = (l-|A|/2), 



A = -1,0,1, 



QA1A2A3 = Aiii'iei + XiKi^i + Asi^aes, 



(9) 
(10) 
(11) 
(12) 

(13) 

(14) 

(15) 

(16) 
(17) 



and {61,62,63} are the unit vectors along the a;, y, and z axes. The Fourier transform of the periodic density has 
max;ima at 



{qi,qi,<l3) = (AiKi, A2iCi, AaisTa), 



(18) 



and each maximum can be labeled by the corresponding set of integers, {A1A2A3}. In units of p(0), the maximum 
of p(q) at q = is 1 (by definition), and the maximum of p(q) at q = (Aii^i, A2-K'2, A3i^3) is approximately 
c(Ai)c(A2)c(A3). The correlation function is 



C{q) = 1 + I /(go) 



^ c(Ai)c(A2)c(A3) 



jl(|q- QaiAsAsI^) 

|q- QA1A2A3I-R 



(19) 



Ai,A2,A3=— 1 

In this case, [C{q) — 1] factorizes, and we can introduce a three-dimensional correlation function C(q) defined by 

C{q)-l=\f{qomC{cD-l]. (20) 

For brevity of notation, we have used the same symbol C for both the three-dimensional and the four-dimensional 
correlation functions in the above equation. The ambiguity of the meaning of C can be easily resolved by context and 
by its argument. 

The three-dimensional correlation function C(q) has maxima values at the same locations as /5(q): at the origin 
q = and at locations qj = ±2Tr\Xj\/dj. When R » d, the overlap of the individual distribution at various maxima 
is small and we have approximately 



C(q)~l+ Yl 



Al ,-^2 ,-^3 



c(Ai)c(A2)c(A3) 



ji(|q- QaiA2A3|^) 
|q- QA1A2A3I-R 



(21) 



The presence of this type of maxima in the three-dimensional correlation function C(q) indicates that prominent 
fluctuations are expected for a spatially periodic granular structure. 
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IV. GENERAL GRANULAR STRUCTURE 



The occurrence of a first-order phase transition in the high temperature environment of a high-energy heavy-ion 
collision will generally lead to granular density distributions that are more general than those discussed in the last 
section. One can consider a general granular density distribution of iV droplets of the type 



N 



(22) 



where pj is the density distribution of the j-th droplet, Xj = (Tj,Rj) is the localized space-time coordinate of the 
center of the j-th droplet, and A is a normalization constant such that the total density p{x) is normalized to unity 
as in Eq. ©. 

We again assume that the source is chaotic in nature. The two-boson correlation function C{q) is then related 
simply to the Fourier transform of the source density. The latter quantity can be easily evaluated and found to be 



p{q) = Jdx e''-p(x) 



N 

where Pj{q) is the Fourier transform of the density distribution of the j-th droplet, 
The correlation function is then 



Ciq) = 1 + 



N 



(23) 



(24) 



(25) 



Thus, the correlation function of a granular source is related to the absolute square of the coherent sum of the Fourier 
transforms of the droplets modulated by the phase e*^ '^^ , which depends on the droplet center coordinate Xj. 

We study the case where the space-time density distribution of a droplet is given by a normalized Gaussian distri- 



bution with standard deviations aj and tj, 



Pjit,r) 



(%/27rcrj)3 VSttTj- 



The Fourier transform of the density of a single droplet is 

p,iq) = J dtdre^**~''i >,(i,r) 
The total density of the granular source is the sum of the density of the droplets. 



e ' e J 



N 



1 



-exp 



(r-R,)^ (t-T,y 



(V27ra,)V2^T, 

The correlation function for the Gaussian granular droplets is then 



2af 



(26) 



(27) 



(28) 



C(g) = 1 + 



N 

AY^e 



(29) 



In g eneral, the droplet lifetime Tj may depend on the droplet size parameter aj (as in the hydrodynamical model 
[33), and (Tj can be different for different droplets. 
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In order to get a clear insight into the most important features of the correlation function of a granular structure, 
we consider the simple case where the density distributions of all droplets are the same so that Uj = Ud and Tj=Td 
for all j. In this simple case, A — l/N and the correlation function can be easily evaluated in terms of the positions 
of the droplet centers, 



This leads to the simple result 

Ciq) 



C{q) = 1 + 



2 2 2 2 



2 2 2 2 



iV2 



N 



N 



N - 



(30) 



(31) 



{Xj — Xk) ~ 2n7r, with n = 1,2,3, .... It has 



Thus, the correlation function C{q) has maxima at q = and qyj 
minima at g^in • [Xj — Xk) ~ (2n — 1)tt. 

We shall consider further the simplified case in which the droplet emission times Tj are the same, then we have 



Ciq) - 1 



2 2 2 2 



iV2 



N 



(32) 



The function [C(g) — 1] factorizes and we can introduce a three-dimensional correlation function C(q) defined by 

C{q)-l = e-^'S'^»[C{q)-l]. (33) 

We have again used the same symbol C for both the three-dimensional and the four-dimensional correlation function 
in the above equation. From Eq. (|30|l . the three-dimensional correlation function C(q) is given by 



C(q) = 1 



7V2 



(34) 



This leads to the simple result 



C(q) = 1 



2 2 



N + 2 



N 

E 

j,k=l,j>k 



cos{q • (Rj - Rfc)} 



(35) 



The correlation function C(q) has a maximum at q = 0, and C(q)|q=o — 2. This maximum at q = is a common 
maximum as it occurs in all correlation functions of two identical bosons. From Eq. H35|l (and if correlation function 
near the maxima and minima does not overlap), the correlation function C(q) will have maxima approximately at 



qmax(i^, n) |2n7r - sin 
and C(q) will have minima approximately at 



{N + 2)aj 



for j,k—l,...N, and j ^ k, 



(36) 



qmin(jfc, ") ^ < (2ri - l)7r - sin 



iN-2)aj 



^5^, for j, fc = 1, ...N, and j ^ k, 

^3k 



(37) 



where n = 1, 2, , 



R 



Ri 



Rfc, and Rjk — |Rj7c|- These maxima and minima correspond to locations where the 



cosine function in Eq. (|35|l has the values close to -f 1 and —1 respectively. 

It is easy to understand how these local maxima and minima of the correlation function of a granular structure 
arise. The interference of the histories of two bosons having a relative momentum q and originating from droplets 
j and k leads to a phase difference of q • (Rj — R^). This phase interference gives rise to the cosine function in Eq. 
(|35|l and a local maximum of the correlation function when the interference is constructive, with a phase difference 
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close to 2n7r. It gives rise to a local minimum when the interference is destructive, with a phase difference close to 
(2n- l)n. 

For each value of n, there are N{N—1) maxima and N{N— 1) minima of the correlation function C(q) of a granular 
structure, one for each permutation of the pairs of droplets. For each maximum located at qmax(j^, n), there is another 
maximum located in the opposite direction at qmax(fci, n) — — qmax(jft, n). Similarly, qmin(fcj, n) = — qminOfc, n). The 
vectors qmax(fcj, n) and qmin(jfc, n) are along the direction of the relative coordinate Rj — R^, and their magnitudes are 
inversely proportional to the magnitudes of the relative coordinates. The magnitude of the maximum of the correlation 
function [C(q) — 1], corresponding to the pair of droplet centers j and fc, is approximately e^'^^'ini^xW'''^) (^]\[ + 2)/N^, 
and the magnitude of the minimum is approximately e~ "''''^"'"^^^ (N — 2)/N'^. 

The maxima and minima of the correlation function C (q) of a granular structure, as given approximately by Eq. 
(|36|l and H37I) . will maintain its distinct characteristics if they are well separated. However, these maxima and minima 
associated with different pairs of droplets may be located in the vicinity of each other. They will merge to give rise 
to a more complicated pattern. When the relative coordinates of many pairs of droplet centers are approximately the 
same (as in the example of the periodic granular structure discussed in Section III), the corresponding maximum of 
the correlation function will be enhanced. The shape of the correlation function will be modified in the presence of 
overlapping maxima and minima, to make the individual maximum and minimum less distinct. It is of interest to 
study some numerical examples to see whether these interesting features of maxima and minima in the correlation 
function manifest themselves. 

If the coordinates of the droplet centers {Rj, j — 1, ...A^} are known, then the correlation function C(q) can 
be evaluated numerically. We consider the center of a droplet ( j-th droplet, say) in many different events to be 
distributed according to a probability distribution P(Rj), 

dP = P{Rj)dRj, (38) 

and assume P(Rj) to be a normalized Gaussian distribution with a standard deviation afj, 

dP = dR,-. (39) 

In our numerical example, we randomly select the localized droplet centers according to the distribution of Eq. I|39|) 
with gr = 4 fm and take ad = 1-5 fm for the standard deviation of the droplet. To avoid overlapping droplets, we 
require the droplets be separated by a distance greater than the sum of their root-mean-square radii, 2^/?iad- After 
the positions of the centers of the droplets are selected, we evaluate the correlation function C(q) with Eq. (P5|l . A 
sample result for iV = 4 droplets is shown in Fig. 1, and another sample result for N = 8 droplets is shown in Fig. 2. 
Figs, (la) and (2a) show the spatial configurations of the droplets; Figs, (lb) and (2b) give the correlation functions 
C(q) as a function of and Qy for = 0, Figs. (Ic) and (2c) for = 0.03 GeV/c, and Figs. (Id) and (2d) for 
qz = 0.06 GeV/c. The momenta q^ and qz in Figs. 1 and 2 are in units of GeV/c. 

One observes from Figs. 1 and 2 that there are prominent fluctuations of the correlation function for a density 
distribution of localized droplets. The inversion symmetry C{qx,qy,qz) = C{~qx,—qy,—qz) is present for qz — 
(Figs, (lb) and (2b)). In addition to the maximum of C(|q|) at q = 0, there are maxima and minima at various 
locations of q. The number of maxima for 8 droplets is greater than the number of maxima for 4 droplets. 

The correlation functions for other configurations of 4 and 8 localized droplets exhibit similarly large fluctuations. 
On the average, the smaller the number of droplets, the greater will be the fluctuation. The magnitude of the 
fluctuation decreases as the number of droplets increases, as can be easily deduced from Eq. The presence of 

this type of maxima and minima of a single-event correlation function at many relative momenta is a signature for a 
granular structure and a first-order QCD phase transition. 

Because of the large fluctuations in a first-order phase transition, one expects that the locations of the centers of 
the droplets will be quite different from event to event. The large differences in the locations of the droplet centers 
in different events lead to large differences in the locations of the maxima and minima and large differences in the 
shapes of correlation functions (except for the maximum at q=0). 

It is interesting to examine the correlation function when we average over many events. The correlation function 
C(q) depends on the coordinates of the droplet centers, and the centers have a distribution f (Rj) in different events. 
The average of the correlation function over the different events, (C(q)), is defined as 

N I N 

(C(q)) = j X{dB.,P{B.i)C{q) / j J]dR,P(Rj). (40) 
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From Eq. H34|l and the Gaussian distribution of the droplet centers (|39|l . we obtain 
and we get 

(C(q)) = 1 + le-'^^'i^ + ^e-('^^+'^H)'^^ (42) 

which is identical to the results of Pratt et at. for a pair of bosons from a distributed source 0. Thus, the results 
of Pratt et al. (il for a distributed source is the same as the results of averaging the correlation function over many 
different events. The average correlation function is now a relatively smooth function of |q|, with only minor fluctu- 
ations even for only two and four droplets. The prominent fluctuations that are inherent in single-event correlation 
functions involving the term cosq • (Rj — R^) in Eq. H35() are not present. The large fluctuations are now greatly 
suppressed by the averaging procedure. 

In order to bring out the salient features of the correlation function, we have assumed in Eq. H31|l that the source 
emission times {Tj} for all droplets are the same. This is a reasonable assumption when the phase transition occurs 
over a short duration of time. The source function in Eq. (|28|l then factorizes in spatial and temporal coordinates. 
The correlation function [C{q) — 1] in Eq. (|33|l also factorizes in qo and q, with C(q) given by Eq. 1)35(1 . 

On the other hand, if the phase transition occurs over a long duration, then the emission times {Tj} can be different 
for different droplets. The source function Eq. H28I) cannot be factorized as a product of spatial and temporal functions. 
Consequently, the correlation function [C{q) — 1] cannot be factorized and is given by the general result of Eq. I|31|l . 
The correlation function C'{q) has maxima at q = and gmax ■ [Xj — X^) 2n7r, with n = 1, 2, 3, .... It has minima 
at qmin • {Xj — Xk) ^ (2n — l)7r. The maxima and minima of the correlation function occur at relative momenta go 
related to the relative emission time coordinates, Tj — Tfc. The signature for the granular structure remains distinct, 
if one can make accurate measurements of the correlation function for different cuts in the relative momenta go- 

There are additional complications when one considers the internal hydrodynamical motion and the collective 
motion of the droplets relative to the center of mass. The hydrodynamics of a QGP droplet shows that the freeze-out 
radial coordinate in a QGP droplet is nearly independent of time in a first-order phase transition, when the initial 
energy density of the droplet is only slightly greater than Ec, the QGP energy density at the critical temperature Tc 
(see Fig. 4(b) of Ref. [32]). As a droplet is presumably formed at temperatures near the critical temperature with 
an energy densities close to ec, the assumption of the factorization of each droplet source as a product of spatial and 
temporal functions is reasonable for the case of a quark- gluon plasma with an initial temperature slightly above Tc- 

For the case of a quark- gluon plasma with a high initial temperature much above Tc, the quark-gluon plasma will 
expand and cool. It will make a phase transition when it cools down to the critical temperature with the subsequent 
formation of granular droplets, if the phase transition is first-order in nature. At the moment of phase transition 
at Tc, the newly- formed droplets will acquire an expansion velocity moving away from the center of mass of the 
system. The factorization of the source function as a product of spatial and temporal functions is not possible, and 
the magnitude of the fluctuations of the correlation function may decrease. We shall investigate these effects on the 
correlation function in our future work. 



V. HOW TO INFER THE DENSITY DISTRIBUTION FROM C{q) 

The correlation function C(g) is related to the Fourier transform of the density, p{q). Information on p{x) is encoded 
in C{q). If the correlation function C{q) has been measured experimentally as a function of its relative momentum 
coordinate g, then a proper Fourier transform of the correlation function will provide pertinent information on the 
density distribution p{x). 

We can obtain direct information on the density distribution p{x) by decoding C(g) in the following way. From the 
correlation function C(g), one calculates [C(g) — 1], and one constructs the Fourier transform of [C(g) — I], 

S{x)= J ^e-«-[C(g)-I]. (43) 

In the discussion of the Fourier transform (or the inversion) of the correlation function, the combination of the two 
terms in [C(g) — I] always comes together. For simplicity, we shall often use the term 'the Fourier transform (or the 
inversion) of the correlation function' to mean 'the Fourier transform (or the inversion) of [C(g) — I]'. 
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The function S{x) is the two-particle 'source function' in the Koonin-Pratt formalismfl^ |27| and the imaging 
method of Brown, Danielewicz and their collaborators |U IstI liM |43 | , for the special case for bosons without final- 
state interactions. Much progress has been made in obtaining a representation of this two-particle source function 
S{x) in terms of spherical harmonics [s^ ISTlli^ l43| . For an irregular density distribution and correlation function as 
one encounters in granular droplets, an expansion of S{x) in terms of spherical harmonics will not be adequate. The 
irregularity of the shape of the correlation function as shown in Figs. 1 and 2 calls for a more general method. The 
best method for inverting a general correlation function without symmetry is to use Cartesian coordinates in a three- 
dimensional Fourier transform. We have developed successfully a general three-dimensional Fast Fourier Transform 
method to invert highly irregular three-dimensional correlation functions. As described in details in Appendix A, 
our three-dimensional FFT method consists of performing a sequence of one-dimensional cosine and sine transforms 
in the three coordinate directions. In each of the cosine or sine transforms, we re-arrange the integral so that the 
limits of the integration go from zero to infinity and the integrand possesses the proper symmetric or antisymmetric 
refiectional symmetry, for cosine or sine transform respectively. We test our numerical three-dimensional FFT method 
by applying it to invert a correlation function of granular droplets for which results for the Fourier inversion can be 
easily obtained analytically. 

In our analysis, we focus our attention on the source density function p{x) itself as it directly gives the space-time 
configuration of the source. From the relation between C{q) and p{x) as given by Eq. we get the integral equation 
for the source density function p{x), 

S{x) = [ dx' p{x') p{x' + x). (44) 



We can prove that S{x) possesses inversion symmetry 

S{x) = S{-x). (45) 

The function S{x) is not the source density but is the folding of the source density with itself. To focus our attention 
on the source density p{x) and to emphasize the property of S{x) as the folding of p(x) with p(x), we can call the 
function S[x) alternatively as 'the folding function' of p{x) in the discussion of the source density p(x), in addition to 
the name of 'the source function' in the Koonin-Pratt formalism and imaging methods '33, l43l |44| . The 

folding function S{x) is real and positive definite. The same folding function S{x) is obtained whether one uses e~'"^'^ 
or its complex conjugate e"^'^ in the Fourier transform expression in Eq. H43|l . 

We can illustrate the application of the folding function S{x) with the example of the chaotic source of N Gaussian 
density droplets, Eq. H28I) . studied in the last section. For such a granular source density p{x), the folding function 
S{x) can be obtained analytically. By carrying out the folding integration using Eq. H44I) . the folding function S{x) 
can be easily found to be 

N 



S(x)^A^ V ' c::p ' (R^' ' R^)]^ [t - (T, - T,) 

(V2^a,fe)3V2^r,, \ 2a% 4^ 



(46) 



where cr|j. = o-^ + and rjj. = rj -I- r|. The function S{x) has maxima at r = (Rj — R^) and t ~ (r^ — r/j), in 
addition to the maxima at x = 0. For simplicity, we again assume that aj — ad and Tj — Td for all j . The function 
S{x) for the granular droplets is simplified to be 

Thus, if the folding function near the maxima does not overlap, the folding function S{x) has maxima at locations 
governed by the relative coordinates of the droplet centers. 

We can consider the case when the droplets are all produced at the same time, as for example when the droplets 
are produced at the moment of phase transition. Then the time Tj can be taken to be the same for all j. The four- 
dimensional folding function S{x) factorizes into a three-dimensional part S{t) and a normalized Gaussian distribution 
in time, 

S{x) = S{r)^^, (48) 
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where we use the same symbol S for the three- and four-dimensional folding function. The ambiguity of the meaning 
of S can be easily resolved by context and by its argument. The three-dimensional folding function S(r) is 

S<.r) . ' V e.p(-t(5l^) , (48) 



If the folding function near the maxima does not overlap, the maxima of the three-dimensional function S{r) are 
located at 

r = Rj-Rfe, j,k = l,2,...N. (50) 

In Eq. H49|) for S'(r), there are N terms with j = fc, and these terms contribute additively to the maxima at r = 0. 
The height at the maxima at r = is therefore N times higher than the maximum with j ^ k located at the the 
relative coordinates Rj — R/j. The occurrence of this type of maxima in S(r), in addition to the maxima at r = 0, 
provides another signature for a granular structure of the source and a first-order phase transition of the quark-gluon 
plasma. 

The density distribution we have considered, with both Tj and Tj separately the same for all j, is a special case of 
those density distributions whose spatial and time distributions can be factorized: p{x) = f(t)p{r), where we use the 
same symbol p for the three-dimensional and the four-dimensional density function. For these factorizable density 
distributions, 

C{q)-l = \f{qomCicO-l], (51) 

and 

C(q)-l = |p(q)P, (52) 

where /(go) and p(q) are the Fourier transforms of f{t) and p{r), respectively. Consequently, the function S{x) also 
factorizes and is given by 

S{x)^S{r) [ ^e-'^^^Vigo)]', (53) 



2tt 



where the three-dimensional function S{r) is equal to 



5(r)= J ^e^'i-[C(q)-l]. (54) 

It is of interest to demonstrate the feasibility and the accuracy of the FFT method by using it to invert a correlation 
function and comparing the inversion result with the analytical result. We use the numerical correlation function 
C(q) obtained in our previous examples in Section IV (results as shown in Fig. 1 and 2) as input, and carry out 
the three-dimensional FFT of the correlation function [C(q) — 1] to obtain ^(r), as given in Eq. 1)54(1 . The input 
correlation functions correspond to those of the localized configurations of Figs, (la) and (2a). 

We show in Figs. 3 and 4 the results of the function S{r) at z = 0, 1.94, 3.87, and 5.81 fm obtained by inverting the 
correlation functions using the FFT method. Fig. 3 gives S{r) for the example of 4 droplets of Fig. 1. Fig. 4 gives 
S{r) for the example of 8 droplets of Fig. 2. One observes that besides the maxima at r = 0, the folding function ^(r) 
has many maxima at r = Rj — R^ where j,k = 1, ...N and j ^ k. A granular structure shows up as having many 
maxima in the Fourier transform of [C(q) — 1], in addition to the maxima at r = 0. For the case with 4 droplets 
in Fig. 3, the maxima of S{r) are quite distinctly exhibited. For the case with 8 droplets, the number of maxima 
increases and many maxima merge. However, some individual maxima remain distinctly separated as in Figs. (4b) 
and (4c). 

To study the shape of the function S{x) in more detail, we consider a cut at the plot of Fig. 3 at x = 0, and plot 
S{x = 0, y, z) as a function of y for different values of z in Fig. 5. The results in Fig. (5a) and (5b) have been obtained 
by using the Fast Fourier Transform method for the example of 4 droplets of Fig. 1. Fig. (5a) gives 5(0, y, z) in linear 
scale and (5b) in logarithmic scale. One sees clearly oscillations of the folding function S{x) due to the maxima at 
various Rj — R^ locations. A signature for granular droplets is the presence of this type of maxima of the Fourier 
transform of [C(q) — 1] at various spatial locations. 

We can assess the accuracy of inverting a numerical correlation function with our Fast Fourier Transform method 
by comparing its results with the exact analytical result as given by Eq. (|49|) . The exact analytical results of S'(r) 
versus y for a; = and different z values are shown in linear scale in Fig. (5c) and in logarithmic scale in Fig. (5d). 
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The results from the FFT method match the exact analytical results with a very high degree of accuracy, including 
the detailed shapes of the oscillations and the values of S{r) down to the low density region of y '--^ ±30 fm where 
S(r) is down by 6 orders of magnitude from its maximum value at r = 0. We have successfully developed an accurate 
three-dimensional FFT method to invert a correlation function [C(r) — 1] to obtain its three-dimensional Fourier 
transform S{r). 

The folding fmiction S{r) will be distorted as noises are introduces into the correlation function. The degree of 
distortion will depend on the magnitude of the noise and it would be of interest to see how well the folding function 
can be determined with noises associated with experimental measurements. The high degree of accuracy in the FFT 
method makes it encouraging to apply it to determine the folding function S{r) for the investigation of the source 
density distribution p(r). 

In order to facilitate the application of the three-dimensional Fourier transform using the experimental single-event 
correlation function [C (q) — 1] (or perhaps functions that fit the experimental correlation function) , we give the detailed 
steps of how the three-dimensional Fourier transform can be evaluated in Appendix A. The computer program to 
carry out the three-dimensional Fast Fourier Transform to obtain S{r) from [C(q) — 1] can also be obtained from the 
authors upon request. Brown et al. have pointed out that in practical applications, when the experimental errors are 
larg e, it is important to include the error uncertainties into the equation for the inversion of the correlation function 
|33l 13^ E3j Q| • Brown [s^l also pointed out that when one applies the FFT transform to experimental correlation 
functions, one should take care to treat the experimental error in the measurement by filtering out the noise, and the 
best method is one in which the errors of the measurement are included in the inversion method. 

It is worth emphasizing that Eq. H44|l . which connects the Fourier transform of the correlation fmiction [C{q) — 1] 
to the density function p{x), is a general result. It can be used to obtain other density distributions, in addition to 
the granular density distribution discussed here. Thus, if the correlation function C{q) is experimentally determined, 
one can first evaluate the Fourier transform of [C{q) — 1], which gives the function S{x). The integral equation Eq. 
(|44|1 can then be used to determine the density distribution p{x) by algebraic methods. In the three-dimensional case, 
one can discretize the integral equation Eq. H44|) as 

N 

S{i,j,k) = AxAyAz ^ p{i + i' ,j + j' ,k + k')p{i' ,j' ,k'), (55) 

i' . j' ,k' — 1 

where x = iAx, y = j Ay, and z = kAz. We consider only the region of S{x) and p{x) inside the box of i, j. A: = 1 to 
N and assume that they are zero outside the box. 

With the determination of >S'(r) from a given experimental correlation function by the FFT method, one can in 
principle solved the above equation to obtain the source density distribution p(r). One can for example solve the 
above equation by iteration. In the first iterative step, one uses a initial guessed density distribution p^^^ (r) as one 
of the two density functions in Eq. (|55|l . The equation is then a linear equation of p(r) and can be easily solved. 
The steps to obtain the solution is described in Appendix B. The solution p(r) can then be substituted into Eq. H55() 
to replace one of the two p functions and to continue the iteration. Clearly, the iterative solution p{r) will be the 
solution of Eq. (|44|l or (|55|l if the iteration converges. The successes of the iterative solution will probably depend 
on a good initial guessed solution. It will be of great interest to test how this iterative procedure may be used to 
find a density distribution p(r) for a given experimental correlation function C(q). It is necessary to investigate how 
one can guarantee positive-definite solutions of p(r) in such procedure. Future development to search for methods to 
solve the integral equation Eq. 1)44(1 [or H55|) ] for p{x) from a given S{x) will be of great interest. 

VI. CONCLUSIONS AND DISCUSSIONS 

Recent experiments at RHIC provide ample evidence for a dense matter produced in high-energy heavy ion collisions. 
Is the produced dense matter the quark-gluon plasma? An unambiguous identification of the produced matter as a 
quark-gluon plasma requires the observation of the phase transition from the new form of matter to known hadronic 
matter. 

The signature for a phase transition depends on the order of the phase transition. Witten and many workers noted 
previously that a granular structure of droplets occurs in a first-order QCD phase transition, and the observation of 
the granular structure can be used as a signature for a first-order QCD phase transition HBT interferometry 

is the best experimental tool to examine the space-time density distribution of the produced matter. It can therefore 
be utilized to study the granular structure that occurs in a first-order phase transition of the plasma. 

In the dynamics following a first-order QCD phase transition, the evolved matter will react chemically and thermally. 
It is not known how much the granular density pattern of the phase transition will remain to make it detectable by HBT 
interferometry. It has been argued in conventional theory that HBT interferometry measures the density distribution 
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of the hadron matter at thermal freeze-out, as the re-scattering of bosons is assumed to lead to a chaotic configuration. 
This traditional assumption is subject to question as it was, however, pointed out recentl y th at t he p ropagation of 
bosons in the re-scattering process should be investigated in a quantum description [s^liol liH l42j . Upon using 
the Glauber theory to describe the scattering process, it was found that HBT interferometry measures the initial 
chaotic density distribution modified by absorption and collective flow The HBT interferometry may be 

sensitive to the density distribution that occurs earlier than the thermal freeze-out configuration. If the initial density 
fluctuation is large, a substantial density fluctuation of the granular pattern may remain to make it detectable by 
HBT interferometry. 

Whatever the theoretical predictions on the possible evolution of the produced matter may be, it is ultimately an 
experimental question whether a granular density distribution that occurs at the moment of the phase transition may 
subsequently render itself detectable by the available experimental tool of HBT interferometry. 

We propose new ways to detect a granular density structure using the single-event HBT interferometry. If it can be 
carried out with sufficient accuracy, the single-event HBT interferometry can reveal the density distribution in each 
single event. It can also reveal large fluctuations in the density distribution from event to event, as is expected in a 
first-order phase transition. 

We carry out our analysis with many examples of granular sources. We found that a granular structure is charac- 
terized by large fluctuations of the single-event correlation function. A single-event correlation function has maxima 
and minima at relative momenta that depend on the relative coordinates of the droplet centers. The presence of 
this type of maxima and minima of the single-event correlation function can be used as the signature for a granular 
structure and the first-order QCD phase transition of the quark-gluon plasma. 

If an experimental single-event correlation function is complete and accurate enough, another very simple method 
to search for a granular structure is to take the Fourier transform of the correlation function. The Fourier transform 
of the correlation function leads to the folding of the source density with itself. This Fourier transform possesses 
maxima at spatial coordinates governed by the relative coordinates of the droplet centers. The occurrence of this type 
of maxima in the Fourier transform, in addition to the maxima at r = 0, is another signature for granular droplets 
and a flrst-order quark-gluon plasma phase transition. 

In the present analysis, we have focused our attention on the single-event HBT interferometry to emphasize the 
maximum fluctuations in the correlation function of a granular structure. It will be of great interest to study a few- 
event HBT analysis both theoretically and experimentally. The few-event HBT analysis will be necessary for practical 
reasons, if the number of boson pairs in a single event is not large enough to provide sufficient statistics. The few-event 
HBT analysis will also be needed to understand the degrees of fluctuation from event to event. One wishes to find 
out whether the fiuctuations in few-event HBT correlation measurements contain a single-event component that is 
beyond statistical fluctuations. The rate of the change of the degrees of fluctuation as the number of events increases 
will provide information on the underlying fluctuations in the single-event and event-to-event HBT interferometry. 
The successful development of the single-event or few-event HBT interferometry in high-energy heavy-ion collisions 
will open up a vast vista for future exploration. 

In the present investigation, we have considered idealized situations in order to bring out the most important 
features of the signature for a granular structure. It will be of great interest to examine in future work how the 
signature discussed here may be affected when some of our simplifying assumptions are modified. The determining 
factor for the occurrence of maxima and minima in the correlation function C (q) is the interference of two histories 
for two bosons to propagate from two source points in different droplets to the detecting points. This interference 
involves a phase difference which depends on the radius vector joining the two droplets. If this underlying factor 
of interference leading to large fluctuations of the correlation function remains important even after modifying our 
simplifying assumptions, then many of the gross features obtained here will not be greatly modifled. 

The fluctuations arising from a granular structure described in the present idealized theoretical investigation will 
be blurred by experimental statistical fluctuations due to the limited number of experimental counts. Whether or 
not the relevant signal can be recovered in the presence of experimental statistical fluctuations remains to be tested. 
Clearly, the smaller the number of droplets, the greater is the signal and the greater will be the probability of its 
observation in the presence of statistical fluctuations. It will be of great interest to carry out a theoretical simulation 
to see what are the minimum droplet size and the largest droplet number a given experimental arrangement may be 
able to detect. 

Much work remains to be done both experimentally and theoretically to investigate this interesting topic on the 
signature of the phase transition of the quark-gluon plasma. It will be of interest to study theoretically effects of 
the collective expansion of the sources and the droplets, effects of fluctuations of the size of the droplets, effects of 
absorption of the bosons on its way to the detector, effects of the momentum dependences of the density distributions, 
and other interesting questions in connection with the signature of the granular structure. 
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APPENDIX A: EVALUATION OF THE THREE-DIMENSIONAL FOURIER TRANSFORM OF THE 

CORRELATION FUNCTION [C(q) - 1)] 



From Eq. (|54(l . we have 

= /(^Wq.r)+^sin(q.r)}[C(q)-l]. (Al) 

The imaginary part in the above integration vanishes as C(— q) — C(q). It is only necessary to evaluate the real part 
of the Fourier transform. The folding function S{r) is therefore 

'S'W^y^^^ cos{qa;X + qyy + Qzz) R{q), (A2) 

where -R(q) = C(q) — 1- Expanding the cosine function, we obtain 

dci . . . 

^^^^3 [cosqxxcosqyy cos qzZ - cosq^xsmqyy smq^z 

— sinq^x cos qyysinqzZ ~ sinq^x sinqyy cos q^z] R{q). (A3) 



In the standard numerical subroutines such as those in the Fast Fourier Transform package of DFFTPACK [6y| , the 
cosine transform of an even function of G(q) is usually approximated by 

G{x)^ d(z cos(gr)G(g) - Ag^cos{(j-l)(fc-l)7r/iV}G(g), (A4) 

where x = (fc — l)Aa:, q — [j ~ 1)A(7, and (Aa:)(Ag) = tt/N . The right-hand side quantities are then calculated by 
the one-dimensional cosine Fast Fourier Transform subroutine of the package. Similarly, the sine transform of an odd 
function is usually approximated by 

G{x)^ sin(gr)G(g) = Ag Vsin{jfc7r/(Af + l)}G(g), (A5) 

with x = kAx and q = jAq and AxAq = it / {N + 1). 

The Fourier integral of the sine and cosine function in Eq. (|A3|I can be cast into the standard form of Eqs. (|A4|I 
and lHU in the EFT package of DFFTPACK by noting that 



dq,^ cos q^x g{qx,qy,qz) = / dq^; cosqxx[g{qx,qy,qz) + g{-qx,qy,qz)], (A6) 

J ~oo Jo 

and similarly 

/OO pOG 
dqx sin q^x g{qx,qy,qz) ^ / dq^ sin q^x [g{qx, qy, qz) - g{-qx,qy,qz)]- (AT) 
-OO Jo 

Using the above results, the four terms inside the square bracket in Eq. IjASp lead to four contributions to S{r), 

Sir) = [Ai(r) - A^ir) - A,iv) - A,iv)] / i27r)\ (A8) 
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where 

Ai{x,y,z)=l dqz cosq^z[Fi{x,y,q:,) + Fi{x,y,-qz)], (A9) 







Fi{x,y,qz)^l dqy cosqyy[Ei{x,qy,qz) + Ei{x,-qy,qz)], (AlO) 
Jo 

and 

Eiix,qy,q,J^ dq^ cosq^x[R{q^,qy,q;,) + R{-qx,qy,q,_)] . (All) 
Jo 

The term A2 that contributes to S{x) can be obtained similarly as 

/"OO 

A2{x,y,z)= / dq^ sinq^z[F2(a;, y, g^) - i^2(a;, -Q^)], (A12) 
^0 



F2ix,y,q^)^ dqy smqyy[E2{x,qy,q^) ~ E2ix,~qy,q^)], (A13) 
Jo 

and 

E2{x,qy,qz)= dqx cos q^x[R{qx,qy,qz) + R{-qx,qy,qz)] ■ (A14) 
Jo 

In a similar way, the term A3 (a;, y, z) is given by 

A3{x,y,z)= / dqz sin g^z[i^3(a;, y, g^) - ^3(2:, ?/, -g^)], (A15) 
Jo 



Fz{x,y,qz)^ i dqy cosqyy[E3{x,qy,qz) + Esix^-qy^q^)], (AI6) 
Jo 

and 

/•OO 

E3{x,qy,qz)^ / dga; sin(gfa;x) [i?(ga;,gs/,9z) - ■ (A17) 
Jo 

Finally, Ai{x,y, z) is given by 

/■OO 

A4(x,?/,z)= / cosqzz[F4{x,y,qz) + F4{x,y,~qz)], (AI8) 

Jo 







FA{x,y,qz)=l dqy sm.qyy[Ei{x,qy,qz) ~ Ei{x,~qy,qz)\, (A19) 

and 

Ei{x,qy,qz) ^ I dqx sinqr^x [R{qx, qy, q^) - R{-qx, qy, qz)] ■ (A20) 





The right-hand sides of the above equations ljA9p - ljA20|) are now in the form of the sine and cosine integrals of Eqs. 
(jA4|) and ljA5p . for which standard FFT subroutines can be applied. With the above relations, the folding function 
S{r) can be easily evaluated using subroutines in standard Fast Fourier Transform packages. 

Incidentally, we have shown how we can obtain the 3-dimensional Fourier transform for a function that is symmetric 
with respect to the inversion of its coordinates, for which the imaginary part of the Fourier transform vanishes. For 
other general functions without such a symmetry, the imaginary Fourier component involving sin(q • r) does not 
vanish. We can expand the function sin.{qxX + qyy + qzz) as in Eq. (|A3() and use techniques similar to those in Eqs. 
^20|l to get the imaginary part of the Fourier transform. 
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APPENDIX B: SOLUTION OF THE DISCRETIZED INTEGRAL EQUATION 

We wish to obtain an iterative solution of the discretized integral equation H55|l for p{v), with a given S{r). As the 
density function p(r) is zero outside the box of {i,i,k — 1 to N}. the summation in Eq. (|55|l can be limited to density 
functions inside the box. Consequently, Eq. H55|) contains fewer and fewer numbers of unknown variables of p as the 
indices i, j, or k of S{i,j,k) increases. We can choose our starting point to be the linear equation containing only 
one variable. The equation can be easily solved. The variables in the subsequent set of linear equation can be solved 
in sequence. Similar procedures can be carried out in two and three dimensions. We shall give the detail procedures 
below for the one-dimensional case to indicate how the iterative solution can be obtained. 

We seek an iterative solution of p{i) satisfying 

N 

5(z)/Ax = 5]p(z + Op("nO, (Bl) 

i' = l 

where p^^"^ («') is either a guessed solution or the solution of the previous iteration. 

Because the density function is zero outside the region of i < 1 and i > N, the set of equations of ljBl|l are 

S{N) = 0, (B2) 



S{N^l)/Ax^p{N)p^°Hl), (B3) 
S{N - 2)/ Ax = p{N - l)p^"\l) + p(7V)p(0)(2), (B4) 
S{N - 3)/Ax = p{N - 2)p(")(l) + p{N - l)p^"\2) + p{N)p^°\3), .... (B5) 

S{l)/Ax = p{l)p("\l) + p(2)p(°)(2) + ... + p{N)p^"Hn). (B6) 

Eq. (|B3|I contains only a single unknown, p{N), and can be solved in terms of the other known quantities. Knowing the 
value of p{N), Eq. (|B4|I can be solved for p{N ~ 1), and so on. In this way, the whole array of p{i) can be determined. 
The above method can be easily generalized to calculate the iterative solution of the 3-dimensional density function 
p{r) from S{r). It will be necessary to normalize the density solution p{r) after each iterative step to ensure that the 
final solution has the proper normalization. 
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APPENDIX C: FIGURES 
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FIG. 1: (a) A sample spatial configuration of four droplets. We use dashed lines to join the droplet centers and the origin to 
show the locations of the droplets. We use vertical dashed lines to indicate the projection of the droplet centers onto the x — y 
plane, (b) Two-boson correlation function C{qx,qy,qz) for — 0, (c) for q^ — 0.030 GeV/c, and (d) q^ = 0.06 GeV/c. The 
quantities qx and qy are in units of GeV/c. 
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(c) g;,=0.03 GeV/c {d) q,=Om GeV/c 



FIG. 2: (a) A sample configuration of eight droplets. We use dashed lines to join the droplet centers and the origin to show 
the locations of the droplets. We use vertical dashed lines to indicate the projection of the droplet centers onto the x — y plane, 
(b) Two-boson correlation function C{qx,qy,qz) as a function of qx and qy at q^ = 0, (c) q^ — 0.03 GeV/c, and (d) qz = 0.06 
GeV/c. The quantities qx and qy are in units of GeV/c. 
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FIG. 3: (a) The folding function S{x,y,z) for 4 droplets in units of fm ^ for the droplet configuration of Fig. 1. Fig. (3a) is 
for 5'(r) at z = 0, (3b) for ^(r) at 2 = 1.94 fm, (3c) for ^(r) at 2 = 3.87 fm, and (3d) for ^(r) at 2 = 5.81 fm. 
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(c) z = 3.87 fm (d) z = 5.81 fm 

FIG. 4: (a) The function S{x,y,z) for 8 droplets in units of fm~^ for the droplet distributions of Fig. 2. Fig. (4a) is for 5(r) 
at 2 = 0, (4b) for ^(r) at z = 1.94 fm, (4c) for ^(r) at z = 3.87 fm, and (4d) for S{r) at z = 5.81 fm. 
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FIG. 5: The function S{x = 0, y, z) versus y for various values of z for the 4 droplets examples of Figs. 1 and 3. The results 
from the FFT Method are given in linear scale in Fig. (5a) and in logarithmic scale in Fig. (5b). The results from the exact 
analytical solution of Eq. H49|l are given in linear scale in Fig. (5c) and in logarithmic scale in Fig. (5d). 



